function res = hydro(GRID_SIZE, L, Z_SIZE, Dt2, tfinal, GAUSS_AMP, GAUSS_SIZE)
  % switches
  HORIZON_ON = false;
  METRIC_ON = false;

  SHOW_PLOT_HYDRO   = true;
  SHOW_PLOT_METRIC  = false;
  SHOW_PLOT_HORIZON = false;
 


  % dirty tricks to make arguments global
  global Dt DX DZ;
  Dt = Dt2;

  TEMP_OFFSET =1.0;

  % First order hydro
  %  eliminating the u-dot-dot terms 
  %  fsolve search method

  % prepare figures for plotting
  close all;

  % load color maps
  load('cm', 'mycmap');
  CMAP_OFFSET = 512;

  if (SHOW_PLOT_HYDRO)
    fighydro  = figure; 
    set(fighydro, 'Colormap', mycmap);
  end
  if (SHOW_PLOT_METRIC)
    figmetric = figure;
    set(figmetric, 'Colormap', mycmap);
  end
  if (SHOW_PLOT_HORIZON)
    fighorizon = figure;
  end


  % compute constants & the r dimension
  DX = L / GRID_SIZE;
  TIME_STEPS = round(tfinal / Dt);
  VECTOR_SUBSAMPLING = 16;
  [gridx, gridy] = meshgrid([1:VECTOR_SUBSAMPLING:GRID_SIZE], [1:VECTOR_SUBSAMPLING:GRID_SIZE]);

  % reserve space for the Gdd computation
  ressize = TIME_STEPS*GRID_SIZE*GRID_SIZE*Z_SIZE*4*4*4;
  %disp(sprintf('Data size   : %.2fMB', ressize/1024/1024));
  %res = zeros(TIME_STEPS, GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
  res = zeros(TIME_STEPS, GRID_SIZE*GRID_SIZE*5);  % integrated Ricci scalar

  % R --- coordinate values
  Z_CUTOFF = 1/(1+TEMP_OFFSET/2);
  DZ = Z_CUTOFF/(Z_SIZE - 1);
  global Z;
  Z = linspace(0, Z_CUTOFF, Z_SIZE);

  % summary
  disp(sprintf('Grid Size   : %d', GRID_SIZE));
  disp(sprintf('Grid Length : %f', L));
  disp(sprintf('Cell Length : %f', DX));
  disp(sprintf('Target Time : %f', tfinal));
  disp(sprintf('Time step   : %f', Dt));
  disp(sprintf('Total steps : %d', TIME_STEPS));
  disp(sprintf('Largest array  (X*X*Z*4*4*4) : %fMB', GRID_SIZE*GRID_SIZE*Z_SIZE*4*4*4/1024/1024));
  pause(1);


  % video stuff
  VIDEO = false;
  EVERY_FRAME=true;
  FPS = 30;
  SECONDS_PER_LIGHTYEAR = 5;
  if (VIDEO)
    vt = VideoWriter('metric.avi', 'Motion JPEG AVI');
    vt.Quality = 100;
    vt.FrameRate = 30;
    open(vt);
  end




  global A;
  A = sqrt(3)*pi/18 - log(3)/2;

  T  = zeros(GRID_SIZE, GRID_SIZE);
  ux = zeros(GRID_SIZE, GRID_SIZE);
  uy = zeros(GRID_SIZE, GRID_SIZE);
  ut =  ones(GRID_SIZE, GRID_SIZE);
  P1 = zeros(GRID_SIZE, GRID_SIZE);
  P2 = zeros(GRID_SIZE, GRID_SIZE);

  
  Tdt  = zeros(GRID_SIZE, GRID_SIZE); 
  uxdt = zeros(GRID_SIZE, GRID_SIZE); 
  uydt = zeros(GRID_SIZE, GRID_SIZE); 
  utdt = zeros(GRID_SIZE, GRID_SIZE);
  P1dt = zeros(GRID_SIZE, GRID_SIZE);
  P2dt = zeros(GRID_SIZE, GRID_SIZE);

  Tdx = zeros(GRID_SIZE, GRID_SIZE); 
  Tdy = zeros(GRID_SIZE, GRID_SIZE); 
  uxdx = zeros(GRID_SIZE, GRID_SIZE); 
  uxdy = zeros(GRID_SIZE, GRID_SIZE); 
  uydx = zeros(GRID_SIZE, GRID_SIZE); 
  uydy = zeros(GRID_SIZE, GRID_SIZE); 
  P1dx = zeros(GRID_SIZE, GRID_SIZE); 
  P1dy = zeros(GRID_SIZE, GRID_SIZE); 
  P2dx = zeros(GRID_SIZE, GRID_SIZE); 
  P2dy = zeros(GRID_SIZE, GRID_SIZE); 

  u = zeros(GRID_SIZE, GRID_SIZE, 3, 4);

  utdtdt = zeros(GRID_SIZE, GRID_SIZE);
  utdtdx = zeros(GRID_SIZE, GRID_SIZE);
  utdtdy = zeros(GRID_SIZE, GRID_SIZE);
  utdxdx = zeros(GRID_SIZE, GRID_SIZE);
  utdxdy = zeros(GRID_SIZE, GRID_SIZE);
  utdydy = zeros(GRID_SIZE, GRID_SIZE);

  uxdtdt = zeros(GRID_SIZE, GRID_SIZE);
  uxdtdx = zeros(GRID_SIZE, GRID_SIZE);
  uxdtdy = zeros(GRID_SIZE, GRID_SIZE);
  uxdxdx = zeros(GRID_SIZE, GRID_SIZE);
  uxdxdy = zeros(GRID_SIZE, GRID_SIZE);
  uxdydy = zeros(GRID_SIZE, GRID_SIZE);

  uydtdt = zeros(GRID_SIZE, GRID_SIZE);
  uydtdx = zeros(GRID_SIZE, GRID_SIZE);
  uydtdy = zeros(GRID_SIZE, GRID_SIZE);
  uydxdx = zeros(GRID_SIZE, GRID_SIZE);
  uydxdy = zeros(GRID_SIZE, GRID_SIZE);
  uydydy = zeros(GRID_SIZE, GRID_SIZE);

  global g gd g_old gi;
  global detg gr;
  global gam gamd gam_old;
  global GUd Gdd;
  g       = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
  g_old   = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
  detg    = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE);
  gr      = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE);
  gam     = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4, 4);
  gam_old = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4, 4);
  GUd     = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);
  Gdd     = zeros(GRID_SIZE, GRID_SIZE, Z_SIZE, 4, 4);





  sumT = 0;
  % initial conditions for vx vy come here
  % vx
  % vy
  % T
  cx = 0.5 * GRID_SIZE + 0.5;
  cy = 0.5 * GRID_SIZE + 0.5;
  gaussian_size = GAUSS_SIZE * L;
  sig2 = (gaussian_size/DX).^2 ;
  for x = 1:GRID_SIZE
    for y = 1:GRID_SIZE
      %T(x, y) =  GAUSS_AMP .* exp(-((x-cx)^2 + (y-cy)^2)/sig2);
      ux(x,y) = ux(x,y) -GAUSS_AMP.*(y-cy).*exp(-((x-cx)^2 + (y-cy)^2)/sig2);
      uy(x,y) = uy(x,y) +GAUSS_AMP.*(x-cx).*exp(-((x-cx)^2 + (y-cy)^2)/sig2);
      %T(x,y) = 0.2*sin(2*pi*x/GRID_SIZE);
    end
  end
  
  T = T + TEMP_OFFSET;



  % ****************************************************************
  %  compute derivative kernels
  % ****************************************************************
  % can upgrade on this -- trefethen pg 6

  global d1 d2;
  d1 = zeros(GRID_SIZE, GRID_SIZE);
  d2 = zeros(GRID_SIZE, GRID_SIZE);
  for x = 1:GRID_SIZE
    for y=1:GRID_SIZE
      s1 = 0;
      s2 = 0;
      for k = -GRID_SIZE/2 : GRID_SIZE/2
        s1 = s1 + 2. * pi * i / GRID_SIZE / GRID_SIZE * k * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
        s2 = s2 + (2. * pi * i / GRID_SIZE * k)^2 / GRID_SIZE * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
      end
      d1(x,y) = real(s1)/DX;
      d2(x,y) = real(s2)/DX/DX;
    end
  end
  h = 2*pi/GRID_SIZE;
  col1 = [0 .5*(-1).^(1:GRID_SIZE-1).*cot((1:GRID_SIZE-1)*h/2)]  * h/DX;
  max(max(d1 - toeplitz(col1, col1([1 GRID_SIZE:-1:2]))));
  col2 = [(-pi^2/(3*h*h)-1/6) .5*(-1).^(0:GRID_SIZE-2).*(csc((1:GRID_SIZE-1)*h/2).^2)]  * h*h/(DX*DX);
  % d2 = toeplitz(col2);
  


  % given T ux, uy.

  % compute ut
  ut = sqrt(1+ux.^2+uy.^2);

  % compute spatial derivatives of initial conditions
  utdx = dx(ut);
  utdy = dy(ut);
  utdxdx = dx(utdx);
  utdxdy = dx(utdy);
  utdydy = dy(utdy);

  uxdx = dx(ux);
  uxdy = dy(ux);
  uxdxdx = dx(uxdx);
  uxdxdy = dx(uxdy);
  uxdydy = dy(uxdy);

  uydx = dx(uy);
  uydy = dy(uy);
  uydxdx = dx(uydx);
  uydxdy = dx(uydy);
  uydydy = dy(uydy);

  Tdx = dx(T);
  Tdy = dy(T);
  Tdxdx = dx(Tdx);
  Tdxdy = dx(Tdy);
  Tdydy = dy(Tdy);

  if false
    % initial conditions for u = 0
    Tdt = 0;
    uxdt = (-1).*T.^(-1).*Tdx;
    uydt = (-1).*T.^(-1).*Tdy;
    utdt = ut.^(-1).*(ux.*uxdt+uy.*uydt);

    Tdtdx = dx(Tdt);
    Tdtdy = dy(Tdt);
    utdtdx = dx(utdt);
    utdtdy = dy(utdt);
    uxdtdx = dx(uxdt);
    uxdtdy = dy(uxdt);
    uydtdx = dx(uydt);
    uydtdy = dy(uydt);

    Tdtdt = (1/2).*(Tdxdx+Tdydy);
    uxdtdt = 0;
    uydtdt = 0;
    utdtdt = 0;

    P1 = (1/4).*(1+A).*T.^(-1).*(2.*Tdx.^2+(-2).*Tdy.^2+T.*((-1).*Tdxdx+Tdydy));
    P2 = (-1/2).*(1+A).*T.^(-1).*(T.*Tdxdy+(-2).*Tdx.*Tdy); 

    P1dx = dx(P1);
    P1dy = dy(P1);
    P2dx = dx(P2);
    P2dy = dy(P2);
  end

  if false

  % compute time derivatives of u and T from 0th order equations of motion.

  % first time derivatives from d_mu T^\mu\nu_0 = 0
  Tdt = (1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).^(-1).*((-1).*Tdx.*ux.*(1+ux.^2)+(-1).*(Tdx.*ux+T.*uxdx).*uy.^2+(-1).*Tdy.*uy.^3+uy.*((-1).*Tdy.*(1+ux.^2)+T.*ux.*(uxdy+uydx))+(-1).*T.*(uxdx+(1+ux.^2).*uydy));
  uxdt = T.^(-1).*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).^(-1).*((-1).*(Tdx+T.*uxdy.*uy).*(2+uy.^2)+(-1).*ux.^2.*(2.*Tdx+T.*uy.*(2.*uxdy+uydx))+(-1).*ux.*(Tdy.*uy+T.*(uxdx+(-1).*uydy))+T.*ux.^3.*((-1).*uxdx+uydy));
  uydt = T.^(-1).*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).^(-1).*((-1).*Tdx.*ux.*uy+(-1).*Tdy.*(2+ux.^2+2.*uy.^2)+(-1).*T.*(ux.^3.*uydx+ux.*(2.*uydx+uy.^2.*(uxdy+2.*uydx))+uy.*(1+uy.^2).*((-1).*uxdx+uydy)));
  utdt = ut.^(-1).*(ux.*uxdt+uy.*uydt);

  % mixed second order
  % the T mixed derivatives are used only in this section
  Tdtdx = dx(Tdt);
  Tdtdy = dy(Tdt);
  utdtdx = dx(utdt);
  utdtdy = dy(utdt);
  uxdtdx = dx(uxdt);
  uxdtdy = dy(uxdt);
  uydtdx = dx(uydt);
  uydtdy = dy(uydt);

  % second order by differentiating the above expressions.
  uxdtdt = (1/2).*T.^(-2).*(ut.^2).^(-3/2).*(2+ux.^2+uy.^2).^(-2).*((-2).*Tdt.*ut.^2.*(2+ux.^2+uy.^2).*((-1).*(Tdx+T.*uxdy.*uy).*(2+uy.^2)+(-1).*ux.^2.*(2.*Tdx+T.*uy.*(2.*uxdy+uydx))+(-1).*ux.*(Tdy.*uy+T.*(uxdx+(-1).*uydy))+T.*ux.^3.*((-1).*uxdx+uydy))+(-4).*T.*ut.^2.*(ux.*uxdt+uy.*uydt).*((-1).*(Tdx+T.*uxdy.*uy).*(2+uy.^2)+(-1).*ux.^2.*(2.*Tdx+T.*uy.*(2.*uxdy+uydx))+(-1).*ux.*(Tdy.*uy+T.*(uxdx+(-1).*uydy))+T.*ux.^3.*((-1).*uxdx+uydy))+(-2).*T.*(2+ux.^2+uy.^2).*(ux.*uxdt+uy.*uydt).*((-1).*(Tdx+T.*uxdy.*uy).*(2+uy.^2)+(-1).*ux.^2.*(2.*Tdx+T.*uy.*(2.*uxdy+uydx))+(-1).*ux.*(Tdy.*uy+T.*(uxdx+(-1).*uydy))+T.*ux.^3.*((-1).*uxdx+uydy))+2.*T.*ut.^2.*(2+ux.^2+uy.^2).*((-2).*uy.*(Tdx+T.*uxdy.*uy).*uydt+(-1).*(2+uy.^2).*(Tdtdx+(T.*uxdtdy+Tdt.*uxdy).*uy+T.*uxdy.*uydt)+T.*ux.^3.*((-1).*uxdtdx+uydtdy)+(-2).*ux.*uxdt.*(2.*Tdx+T.*uy.*(2.*uxdy+uydx))+(-1).*ux.^2.*(2.*Tdtdx+T.*uy.*(2.*uxdtdy+uydtdx)+Tdt.*uy.*(2.*uxdy+uydx)+T.*uydt.*(2.*uxdy+uydx))+(-1).*uxdt.*(Tdy.*uy+T.*(uxdx+(-1).*uydy))+(-1).*ux.*(Tdtdy.*uy+Tdy.*uydt+T.*(uxdtdx+(-1).*uydtdy)+Tdt.*(uxdx+(-1).*uydy))+Tdt.*ux.^3.*((-1).*uxdx+uydy)+3.*T.*ux.^2.*uxdt.*((-1).*uxdx+uydy)));
  uydtdt = (1/2).*T.^(-2).*(ut.^2).^(-3/2).*(2+ux.^2+uy.^2).^(-2).*((-2).*Tdt.*ut.^2.*(2+ux.^2+uy.^2).*((-1).*Tdx.*ux.*uy+(-1).*Tdy.*(2+ux.^2+2.*uy.^2)+(-1).*T.*(ux.^3.*uydx+ux.*(2.*uydx+uy.^2.*(uxdy+2.*uydx))+uy.*(1+uy.^2).*((-1).*uxdx+uydy)))+(-4).*T.*ut.^2.*(ux.*uxdt+uy.*uydt).*((-1).*Tdx.*ux.*uy+(-1).*Tdy.*(2+ux.^2+2.*uy.^2)+(-1).*T.*(ux.^3.*uydx+ux.*(2.*uydx+uy.^2.*(uxdy+2.*uydx))+uy.*(1+uy.^2).*((-1).*uxdx+uydy)))+(-2).*T.*(2+ux.^2+uy.^2).*(ux.*uxdt+uy.*uydt).*((-1).*Tdx.*ux.*uy+(-1).*Tdy.*(2+ux.^2+2.*uy.^2)+(-1).*T.*(ux.^3.*uydx+ux.*(2.*uydx+uy.^2.*(uxdy+2.*uydx))+uy.*(1+uy.^2).*((-1).*uxdx+uydy)))+2.*T.*ut.^2.*(2+ux.^2+uy.^2).*((-1).*Tdtdx.*ux.*uy+(-1).*Tdx.*uxdt.*uy+(-1).*Tdtdy.*(2+ux.^2+2.*uy.^2)+(-1).*Tdx.*ux.*uydt+(-2).*Tdy.*(ux.*uxdt+2.*uy.*uydt)+(-1).*Tdt.*(ux.^3.*uydx+ux.*(2.*uydx+uy.^2.*(uxdy+2.*uydx))+uy.*(1+uy.^2).*((-1).*uxdx+uydy))+(-1).*T.*(ux.^3.*uydtdx+uy.*(1+uy.^2).*((-1).*uxdtdx+uydtdy)+3.*ux.^2.*uxdt.*uydx+uxdt.*(2.*uydx+uy.^2.*(uxdy+2.*uydx))+ux.*(2.*uydtdx+uy.^2.*(uxdtdy+2.*uydtdx)+2.*uy.*uydt.*(uxdy+2.*uydx))+2.*uy.^2.*uydt.*((-1).*uxdx+uydy)+(1+uy.^2).*uydt.*((-1).*uxdx+uydy))));

  utdtdt = ut.^(-3).*(uxdt.^2.*(1+uy.^2)+(-2).*ux.*uxdt.*uy.*uydt+(1+ux.^2).*uydt.^2+ut.^2.*(ux.*uxdtdt+uy.*uydtdt));

  P1 = (1/8).*T.*(1+ux.^2+uy.^2).^(-3/2).*(2+ux.^2+uy.^2).^(-1).*(4.*uxdtdx+4.*A.*uxdtdx+(-4).*(2+A).*ux.^4.*uxdt.*uxdy.*uy.^3+4.*(2+A).*uxdt.*uxdy.*uy.^7+2.*(2+A).*(uxdt.^2+uxdy.^2).*uy.^6.*(1+ux.^2+uy.^2).^(1/2)+2.*(1+A).*ux.^3.*uy.*(7.*uxdtdy+(-9).*uydtdx)+2.*(1+A).*ux.*uy.^3.*(9.*uxdtdy+(-7).*uydtdx)+2.*(1+A).*ux.^5.*uy.*(3.*uxdtdy+(-5).*uydtdx)+2.*(1+A).*ux.*uy.^5.*(5.*uxdtdy+(-3).*uydtdx)+8.*(1+A).*ux.*uy.*(uxdtdy+(-1).*uydtdx)+16.*(1+A).*ux.^3.*uy.^3.*(uxdtdy+(-1).*uydtdx)+(-4).*uydtdy+(-4).*A.*uydtdy+(-2).*(1+A).*ux.^6.*((-2).*uxdtdx+uydtdy)+(-26).*(1+A).*ux.^2.*uy.^2.*((-1).*uxdtdx+uydtdy)+(-2).*(1+A).*uy.^6.*((-1).*uxdtdx+2.*uydtdy)+(-2).*(1+A).*ux.^4.*((-7).*uxdtdx+4.*uydtdy)+(-2).*(1+A).*ux.^2.*((-7).*uxdtdx+5.*uydtdy)+(-2).*(1+A).*ux.^4.*uy.^2.*((-7).*uxdtdx+6.*uydtdy)+(-2).*(1+A).*ux.^2.*uy.^4.*((-6).*uxdtdx+7.*uydtdy)+(-2).*(1+A).*uy.^2.*((-5).*uxdtdx+7.*uydtdy)+(-2).*(1+A).*uy.^4.*((-4).*uxdtdx+7.*uydtdy)+(-4).*(2+A).*ux.^7.*uydt.*uydx+4.*(2+A).*ux.^3.*uy.^4.*uydt.*uydx+(-4).*ux.*(T.*uxdt+(-2).*(1+A).*uxdt.*uxdx+(3+A).*uydt.*uydx)+(-2).*ux.^5.*(T.*uxdt+(-1).*(1+A).*uxdt.*uxdx+(15+7.*A).*uydt.*uydx)+ux.^3.*((-6).*T.*uxdt+8.*(1+A).*uxdt.*uxdx+(-2).*(17+7.*A).*uydt.*uydx)+(-2).*(2+A).*ux.^6.*(1+ux.^2+uy.^2).^(1/2).*(uydt.^2+uydx.^2)+2.*(1+A).*ux.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdtdt+3.*uxdxdx+2.*uxdydy+(-5).*uydxdy)+2.*(1+A).*ux.*uy.^4.*(1+ux.^2+uy.^2).^(1/2).*(3.*uxdtdt+uxdxdx+2.*uxdydy+(-3).*uydxdy)+4.*(1+A).*ux.^3.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdtdt+uxdxdx+uxdydy+(-2).*uydxdy)+(-4).*(1+A).*ux.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdtdt+(-1).*uxdxdx+uydxdy)+(-6).*(1+A).*ux.^3.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdtdt+(-1).*uxdxdx+uydxdy)+(-2).*(1+A).*ux.^5.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdtdt+(-1).*uxdxdx+uydxdy)+4.*(2+A).*ux.^6.*uy.*(uxdt.*uydx+uydt.*(uxdx+(-1).*uydy))+(-4).*(2+A).*ux.^2.*uy.^5.*(uxdt.*uydx+uydt.*(uxdx+(-1).*uydy))+4.*(2+A).*ux.^3.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdy+uydx).*((-1).*uxdx+uydy)+4.*(2+A).*ux.^5.*uy.^2.*((-1).*uxdt.*uxdx+uxdy.*uydt+uxdt.*uydy)+(-4).*(2+A).*ux.*uy.^6.*((-1).*uxdt.*uxdx+uxdy.*uydt+uxdt.*uydy)+(-2).*ux.*uy.^2.*(5.*T.*uxdt+(-13).*uxdt.*uxdx+(-11).*A.*uxdt.*uxdx+uxdy.*uydt+(-1).*A.*uxdy.*uydt+6.*uydt.*uydx+2.*A.*uydt.*uydx+(-1).*((-1)+A).*uxdt.*uydy)+ux.^3.*uy.^2.*((-8).*T.*uxdt+8.*uxdt.*uxdx+8.*A.*uxdt.*uxdx+8.*uxdy.*uydt+8.*A.*uxdy.*uydt+(-12).*uydt.*uydx+(-4).*A.*uydt.*uydx+8.*(1+A).*uxdt.*uydy)+(-2).*ux.*uy.^4.*(3.*T.*uxdt+(-13).*uxdt.*uxdx+(-9).*A.*uxdt.*uxdx+6.*uxdy.*uydt+2.*A.*uxdy.*uydt+uydt.*uydx+A.*uydt.*uydx+2.*(3+A).*uxdt.*uydy)+(-4).*(2+A).*ux.^5.*uy.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdt.*uydt+(-1).*uxdx.*uydx+uydx.*uydy)+(-2).*(2+A).*ux.^2.*uy.^4.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdx.^2+uxdy.^2+(-1).*uydt.^2+2.*uxdy.*uydx+2.*uxdx.*uydy+(-1).*uydy.^2)+(-2).*(2+A).*ux.^4.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.^2+uxdx.^2+(-2).*uxdy.*uydx+(-1).*uydx.^2+(-2).*uxdx.*uydy+uydy.^2)+uy.^4.*(1+ux.^2+uy.^2).^(1/2).*(14.*uxdt.^2+8.*A.*uxdt.^2+uxdx.^2+A.*uxdx.^2+4.*(3+A).*uxdy.^2+(-1).*uydt.^2+(-1).*A.*uydt.^2+(-1).*uydy.^2+(-1).*A.*uydy.^2+2.*T.*((-1).*uxdx+uydy))+2.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdt.^2+2.*A.*uxdt.^2+uxdx.^2+A.*uxdx.^2+(-2).*A.*uxdy.^2+(-2).*uydt.^2+(-2).*A.*uydt.^2+2.*A.*uydx.^2+(-1).*uydy.^2+(-1).*A.*uydy.^2+2.*T.*((-1).*uxdx+uydy))+ux.^4.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.^2+A.*uxdt.^2+uxdx.^2+A.*uxdx.^2+(-14).*uydt.^2+(-8).*A.*uydt.^2+(-12).*uydx.^2+(-4).*A.*uydx.^2+(-1).*(1+A).*uydy.^2+2.*T.*((-1).*uxdx+uydy))+ux.^2.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.^2+3.*A.*uxdt.^2+6.*uxdx.^2+6.*A.*uxdx.^2+(-5).*uydt.^2+(-3).*A.*uydt.^2+(-6).*(1+A).*uydy.^2+4.*T.*((-1).*uxdx+uydy))+uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(14.*uxdt.^2+10.*A.*uxdt.^2+3.*uxdx.^2+3.*A.*uxdx.^2+(-2).*((-4)+A).*uxdy.^2+(-4).*uydt.^2+(-4).*A.*uydt.^2+2.*A.*uydx.^2+(-3).*uydy.^2+(-3).*A.*uydy.^2+6.*T.*((-1).*uxdx+uydy))+ux.^2.*(1+ux.^2+uy.^2).^(1/2).*(4.*uxdt.^2+4.*A.*uxdt.^2+3.*uxdx.^2+3.*A.*uxdx.^2+(-2).*A.*uxdy.^2+(-14).*uydt.^2+(-10).*A.*uydt.^2+(-8).*uydx.^2+2.*A.*uydx.^2+(-3).*uydy.^2+(-3).*A.*uydy.^2+6.*T.*((-1).*uxdx+uydy))+(-4).*(2+A).*ux.*uy.^5.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.*uydt+uxdy.*((-1).*uxdx+uydy))+2.*ux.^3.*uy.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.*uydt+3.*A.*uxdt.*uydt+(-2).*T.*(uxdy+(-1).*uydx)+3.*uxdx.*uydx+(-1).*A.*uxdx.*uydx+(-9).*uydx.*uydy+(-5).*A.*uydx.*uydy+(1+A).*uxdy.*(uxdx+uydy))+(-2).*ux.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.*uydt+3.*A.*uxdt.*uydt+2.*T.*(uxdy+(-1).*uydx)+uxdx.*uydx+A.*uxdx.*uydx+uydx.*uydy+A.*uydx.*uydy+(-1).*uxdy.*((9+5.*A).*uxdx+((-3)+A).*uydy))+uy.^3.*(2.*(17+7.*A).*uxdt.*uxdy+2.*uydt.*(3.*T+(-4).*(1+A).*uydy))+4.*uy.*((3+A).*uxdt.*uxdy+uydt.*(T+(-2).*(1+A).*uydy))+2.*uy.^5.*((15+7.*A).*uxdt.*uxdy+uydt.*(T+(-1).*(1+A).*uydy))+(-2).*ux.*uy.*(1+ux.^2+uy.^2).^(1/2).*(2.*T.*(uxdy+(-1).*uydx)+(-1).*uxdy.*((5+3.*A).*uxdx+(1+3.*A).*uydy)+uydx.*((1+3.*A).*uxdx+(5+3.*A).*uydy))+2.*ux.^4.*uy.*((1+A).*uxdt.*uxdy+2.*(3+A).*uxdt.*uydx+uydt.*(3.*T+2.*(3+A).*uxdx+(-1).*(13+9.*A).*uydy))+2.*ux.^2.*uy.*(2.*(3+A).*uxdt.*uxdy+(-1).*((-1)+A).*uxdt.*uydx+uydt.*(5.*T+(-1).*((-1)+A).*uxdx+(-1).*(13+11.*A).*uydy))+4.*ux.^2.*uy.^3.*((3+A).*uxdt.*uxdy+(-2).*((1+A).*uxdt.*uydx+uydt.*((-1).*T+(1+A).*(uxdx+uydy))))+(-4).*(1+A).*uy.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdxdy+uydtdt+uydydy)+(-6).*(1+A).*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdxdy+uydtdt+uydydy)+(-2).*(1+A).*uy.^5.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdxdy+uydtdt+uydydy)+(-4).*(1+A).*ux.^2.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*((-2).*uxdxdy+2.*uydtdt+uydxdx+uydydy)+(-2).*(1+A).*ux.^4.*uy.*(1+ux.^2+uy.^2).^(1/2).*((-3).*uxdxdy+3.*uydtdt+2.*uydxdx+uydydy)+(-2).*(1+A).*ux.^2.*uy.*(1+ux.^2+uy.^2).^(1/2).*((-5).*uxdxdy+5.*uydtdt+2.*uydxdx+3.*uydydy));

  P2 = (1/4).*T.*(1+ux.^2+uy.^2).^(-3/2).*(2+ux.^2+uy.^2).^(-1).*(2.*uxdtdy+2.*A.*uxdtdy+2.*(1+A).*ux.*uxdtdx.*uy.^5+4.*(1+A).*uxdtdy.*uy.^6+(-4).*(2+A).*ux.*uxdt.*uxdy.*uy.^6+2.*(1+A).*ux.*uxdxdy.*uy.^4.*(1+ux.^2+uy.^2).^(1/2)+(-2).*(2+A).*ux.*(uxdt.^2+uxdy.^2).*uy.^5.*(1+ux.^2+uy.^2).^(1/2)+2.*(1+A).*(uxdtdt+uxdydy).*uy.^5.*(1+ux.^2+uy.^2).^(1/2)+2.*uydtdx+2.*A.*uydtdx+4.*(1+A).*ux.^6.*uydtdx+10.*(1+A).*ux.^2.*uy.^2.*(uxdtdy+uydtdx)+4.*(1+A).*uy.^2.*(2.*uxdtdy+uydtdx)+2.*(1+A).*ux.^2.*uy.^4.*(3.*uxdtdy+uydtdx)+2.*(1+A).*uy.^4.*(5.*uxdtdy+uydtdx)+4.*(1+A).*ux.^2.*(uxdtdy+2.*uydtdx)+2.*(1+A).*ux.^4.*uy.^2.*(uxdtdy+3.*uydtdx)+2.*(1+A).*ux.^4.*(uxdtdy+5.*uydtdx)+2.*(1+A).*ux.^5.*uy.*uydtdy+2.*(1+A).*ux.*uy.*(uxdtdx+uydtdy)+2.*(1+A).*ux.^3.*uy.^3.*(uxdtdx+uydtdy)+2.*(1+A).*ux.*uy.^3.*(2.*uxdtdx+uydtdy)+2.*(1+A).*ux.^3.*uy.*(uxdtdx+2.*uydtdy)+(-4).*(2+A).*ux.^6.*uy.*uydt.*uydx+(-2).*(2+A).*ux.^5.*uy.*(1+ux.^2+uy.^2).^(1/2).*(uydt.^2+uydx.^2)+2.*(1+A).*ux.^5.*(1+ux.^2+uy.^2).^(1/2).*(uydtdt+uydxdx)+2.*(1+A).*ux.*(1+ux.^2+uy.^2).^(1/2).*(uxdxdy+uydtdt+uydxdx)+2.*(1+A).*ux.^3.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(uxdxdy+uydtdt+uydxdx)+2.*(1+A).*ux.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdxdy+uydtdt+uydxdx)+2.*(1+A).*ux.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdxdy+2.*(uydtdt+uydxdx))+2.*(1+A).*ux.^4.*uy.*(1+ux.^2+uy.^2).^(1/2).*uydxdy+2.*(1+A).*uy.*(1+ux.^2+uy.^2).^(1/2).*(uxdtdt+uxdydy+uydxdy)+2.*(1+A).*ux.^2.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdtdt+uxdydy+uydxdy)+2.*(1+A).*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdtdt+2.*uxdydy+uydxdy)+2.*(1+A).*ux.^2.*uy.*(1+ux.^2+uy.^2).^(1/2).*(uxdtdt+uxdydy+2.*uydxdy)+4.*(2+A).*ux.^5.*uy.^2.*(uxdt.*uydx+uydt.*(uxdx+(-1).*uydy))+4.*(2+A).*ux.^2.*uy.^5.*((-1).*uxdt.*uxdx+uxdy.*uydt+uxdt.*uydy)+4.*(2+A).*ux.^4.*uy.^3.*((-1).*uxdt.*uxdx+uydt.*(uxdy+(-1).*uydx)+uxdt.*uydy)+ux.^4.*uy.*((-4).*(2+A).*uxdt.*uxdx+uydt.*(2.*(3+A).*uxdy+(-1).*(19+9.*A).*uydx)+2.*(3+A).*uxdt.*uydy)+uy.^5.*((-2).*T.*uxdt+(-1).*uxdt.*uxdx+A.*uxdt.*uxdx+5.*uxdy.*uydt+3.*A.*uxdy.*uydt+(5+3.*A).*uxdt.*uydy)+uy.*((-2).*T.*uxdt+(-1).*uxdt.*uxdx+A.*uxdt.*uxdx+5.*uxdy.*uydt+3.*A.*uxdy.*uydt+(-1).*uydt.*uydx+A.*uydt.*uydx+(5+3.*A).*uxdt.*uydy)+uy.^3.*((-4).*T.*uxdt+(-2).*uxdt.*uxdx+2.*A.*uxdt.*uxdx+10.*uxdy.*uydt+6.*A.*uxdy.*uydt+(-1).*uydt.*uydx+A.*uydt.*uydx+2.*(5+3.*A).*uxdt.*uydy)+ux.^2.*uy.*((-2).*T.*uxdt+(-11).*uxdt.*uxdx+(-5).*A.*uxdt.*uxdx+11.*uxdy.*uydt+5.*A.*uxdy.*uydt+(-12).*uydt.*uydx+(-4).*A.*uydt.*uydx+(11+5.*A).*uxdt.*uydy)+ux.^2.*uy.^3.*((-2).*T.*uxdt+(-19).*uxdt.*uxdx+(-9).*A.*uxdt.*uxdx+21.*uxdy.*uydt+11.*A.*uxdy.*uydt+(-9).*uydt.*uydx+(-3).*A.*uydt.*uydx+(21+11.*A).*uxdt.*uydy)+(-4).*(2+A).*ux.^4.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdt.*uydt+(-1).*uxdx.*uydx+uydx.*uydy)+ux.^4.*(1+ux.^2+uy.^2).^(1/2).*((-2).*T.*uydx+(5+3.*A).*(uxdt.*uydt+uxdx.*uydx)+((-1)+A).*uydx.*uydy)+(-2).*(2+A).*ux.^3.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.^2+uxdx.^2+uydt.^2+(-2).*uxdy.*uydx+(-2).*uxdx.*uydy+uydy.^2)+(-1).*ux.*uy.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.^2+3.*A.*uxdt.^2+2.*uxdx.^2+2.*uxdy.^2+5.*uydt.^2+3.*A.*uydt.^2+(-4).*uxdy.*uydx+2.*uydx.^2+(-4).*uxdx.*uydy+2.*uydy.^2)+(-1).*ux.^3.*uy.*(1+ux.^2+uy.^2).^(1/2).*(4.*uxdt.^2+2.*A.*uxdt.^2+2.*(2+A).*uxdx.^2+9.*uydt.^2+5.*A.*uydt.^2+(-6).*uxdy.*uydx+(-2).*A.*uxdy.*uydx+6.*uydx.^2+2.*A.*uydx.^2+(-2).*(3+A).*uxdx.*uydy+2.*uydy.^2)+(-1).*ux.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(9.*uxdt.^2+5.*A.*uxdt.^2+2.*uxdx.^2+2.*(3+A).*uxdy.^2+4.*uydt.^2+2.*A.*uydt.^2+(-2).*(3+A).*uxdy.*uydx+(-2).*(3+A).*uxdx.*uydy+2.*(2+A).*uydy.^2)+4.*(2+A).*ux.^2.*uy.^4.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.*uydt+uxdy.*((-1).*uxdx+uydy))+ux.^2.*(1+ux.^2+uy.^2).^(1/2).*(9.*uxdt.*uydt+7.*A.*uxdt.*uydt+6.*uxdx.*uydx+2.*A.*uxdx.*uydx+(-2).*T.*(uxdy+2.*uydx)+4.*A.*uydx.*uydy+(-1).*((-1)+A).*uxdy.*((-1).*uxdx+uydy))+ux.^2.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(20.*uxdt.*uydt+12.*A.*uxdt.*uydt+9.*uxdx.*uydx+3.*A.*uxdx.*uydx+(-2).*T.*(uxdy+uydx)+(-9).*uydx.*uydy+(-3).*A.*uydx.*uydy+3.*(3+A).*uxdy.*((-1).*uxdx+uydy))+(-4).*(2+A).*ux.^3.*uy.^4.*(uxdt.*uxdy+(-1).*uxdt.*uydx+uydt.*((-1).*uxdx+uydy))+(1+ux.^2+uy.^2).^(1/2).*(4.*uxdt.*uydt+4.*A.*uxdt.*uydt+uxdx.*uydx+(-1).*A.*uxdx.*uydx+(-2).*T.*(uxdy+uydx)+uydx.*uydy+3.*A.*uydx.*uydy+uxdy.*((1+3.*A).*uxdx+(-1).*((-1)+A).*uydy))+ux.^5.*((5+3.*A).*uxdt.*uydx+uydt.*((-2).*T+(5+3.*A).*uxdx+((-1)+A).*uydy))+ux.*(((-1)+A).*uxdt.*uxdy+(5+3.*A).*uxdt.*uydx+uydt.*((-2).*T+(5+3.*A).*uxdx+((-1)+A).*uydy))+uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(9.*uxdt.*uydt+7.*A.*uxdt.*uydt+uxdx.*uydx+(-1).*A.*uxdx.*uydx+(-2).*T.*(2.*uxdy+uydx)+(-1).*uydx.*uydy+A.*uydx.*uydy+2.*uxdy.*(2.*A.*uxdx+(3+A).*uydy))+uy.^4.*(1+ux.^2+uy.^2).^(1/2).*((-2).*T.*uxdy+(5+3.*A).*uxdt.*uydt+uxdy.*(((-1)+A).*uxdx+(5+3.*A).*uydy))+ux.^3.*uy.^2.*((-3).*(3+A).*uxdt.*uxdy+(21+11.*A).*uxdt.*uydx+uydt.*((-2).*T+(21+11.*A).*uxdx+(-1).*(19+9.*A).*uydy))+ux.*uy.^2.*((-4).*(3+A).*uxdt.*uxdy+(11+5.*A).*uxdt.*uydx+uydt.*((-2).*T+(-1).*(11+5.*A).*((-1).*uxdx+uydy)))+ux.^3.*(((-1)+A).*uxdt.*uxdy+2.*((5+3.*A).*uxdt.*uydx+uydt.*((-2).*T+(5+3.*A).*uxdx+((-1)+A).*uydy)))+ux.*uy.^4.*((-1).*(19+9.*A).*uxdt.*uxdy+2.*((3+A).*uxdt.*uydx+uydt.*((3+A).*uxdx+(-2).*(2+A).*uydy))));

  P1dx = dx(P1);
  P1dy = dy(P1);
  P2dx = dx(P2);
  P2dy = dy(P2);

  end
  % compute the values of P1 and P2. woooo



  


  % initial conditions for u = 0
  % P1 = (-1/4).*T.^(-1).*((-1).*Tdx.^2+(-1).*A.*Tdx.^2+Tdy.^2+A.*Tdy.^2);
  % P2 = (-1/2).*T.^(-1).*((-1).*Tdx.*Tdy+(-1).*A.*Tdx.*Tdy);

  % P2 = (1/4).*(1+A).*(dy(Tdy)+dx(Tdx));
  % P1 = (1/4).*(1+A).*(dy(Tdy)+dx(Tdx));


    t = 1;
    if (SHOW_PLOT_HYDRO)
      phydro();
      pause(0.003);
    end


  if (METRIC_ON)    
    ut = sqrt(1+ux.^2+uy.^2);
    utdt = zeros(GRID_SIZE, GRID_SIZE);



    Tdx = dx(T);
    Tdy = dy(T);
    utdx = dx(ut);
    utdy = dy(ut);
    uxdx = dx(ux);
    uxdy = dy(ux);
    uydx = dx(uy);
    uydy = dy(uy);
      % compute all second order derivatives

      % back differencing for dtdt
    utdtdt = zeros(GRID_SIZE, GRID_SIZE);
    uxdtdt = zeros(GRID_SIZE, GRID_SIZE);
    uydtdt = zeros(GRID_SIZE, GRID_SIZE);

      % spectral for all others
      utdtdx = dx(utdt);
      uxdtdx = dx(uxdt);
      uydtdx = dx(uydt);

      utdtdy = dy(utdt);
      uxdtdy = dy(uxdt);
      uydtdy = dy(uydt);

      utdxdx = dx(utdx);
      uxdxdx = dx(uxdx);
      uydxdx = dx(uydx);

      utdxdy = dy(utdx);
      uxdxdy = dy(uxdx);
      uydxdy = dy(uydx);

      utdydy = dy(utdy);
      uxdydy = dy(uxdy);
      uydydy = dy(uydy);
       
      % call the metric code
      metricz(T, ut, ux, uy, utdt, uxdt, uydt, utdx, utdy, uxdx, uxdy, uydx, uydy, ...
          utdtdt, utdtdx, utdtdy, utdxdx, utdxdy, utdydy, ...
          uxdtdt, uxdtdx, uxdtdy, uxdxdx, uxdxdy, uxdydy, ...
          uydtdt, uydtdx, uydtdy, uydxdx, uydxdy, uydydy);
  end
  


  for t = 1:TIME_STEPS
    tic;

    % ids of the past values of u
    uid = rem(t+3:-1:t,4) + 1;

    % solving non-linear equation of type:
    % L(u, T, u-dot, T-dot) + R(u, T) == 0  -- 3 components
    % each of L, R has 3 components 

    [Tdt uxdt uydt P1dt P2dt] = h2eom(T, ux, uy, P1, P2, Tdx, Tdy, uxdx, uxdy, uydx, uydy, P1dx, P1dy, P2dx, P2dy);

    % single space derivatives of T, ux, uy, P1, P2 have been computed after previous update
    %  catch up with ut only:
    ut = sqrt(1+ux.^2+uy.^2);
    utdt = (ux.*uxdt + uy.*uydt)./ut;
    utdx = dx(ut);
    utdy = dy(ut);


    % save current values of u for computing of derivatives
    % TODO : this is the missing code
    u(:,:,1,uid(1)) = ut;
    u(:,:,2,uid(1)) = ux;
    u(:,:,3,uid(1)) = uy;
    % force the time derivatives to zero for the very first time step
    if t == 1
      for k=1:4
        u(:,:,1,k) = ut;
        u(:,:,2,k) = ux;
        u(:,:,3,k) = uy;
      end
    end
    %utdt = (u(:,:,1,uid(1)) - u(:,:,1,uid(2)))./Dt;

    % back differencing for dtdt
    % compute second order derivatives 
    utdtdt = (2.*u(:,:,1,uid(1)) - 5.*u(:,:,1,uid(2)) + 4.*u(:,:,1,uid(3)) - 1.*u(:,:,1,uid(4)))./Dt./Dt;
    uxdtdt = (2.*u(:,:,2,uid(1)) - 5.*u(:,:,2,uid(2)) + 4.*u(:,:,2,uid(3)) - 1.*u(:,:,2,uid(4)))./Dt./Dt;
    uydtdt = (2.*u(:,:,3,uid(1)) - 5.*u(:,:,3,uid(2)) + 4.*u(:,:,3,uid(3)) - 1.*u(:,:,3,uid(4)))./Dt./Dt;
    
    % compute all spatial and mixed second order derivatives
    utdtdx = dx(utdt);
    uxdtdx = dx(uxdt);
    uydtdx = dx(uydt);

    utdtdy = dy(utdt);
    uxdtdy = dy(uxdt);
    uydtdy = dy(uydt);

    utdxdx = dx(utdx);
    uxdxdx = dx(uxdx);
    uydxdx = dx(uydx);

    utdxdy = dy(utdx);
    uxdxdy = dy(uxdx);
    uydxdy = dy(uydx);

    utdydy = dy(utdy);
    uxdydy = dy(uxdy);
    uydydy = dy(uydy);

    % check the constitutive relations

    PP1 = (1/8).*T.*(1+ux.^2+uy.^2).^(-3/2).*(2+ux.^2+uy.^2).^(-1).*(4.*uxdtdx+4.*A.*uxdtdx+(-4).*(2+A).*ux.^4.*uxdt.*uxdy.*uy.^3+4.*(2+A).*uxdt.*uxdy.*uy.^7+2.*(2+A).*(uxdt.^2+uxdy.^2).*uy.^6.*(1+ux.^2+uy.^2).^(1/2)+2.*(1+A).*ux.^3.*uy.*(7.*uxdtdy+(-9).*uydtdx)+2.*(1+A).*ux.*uy.^3.*(9.*uxdtdy+(-7).*uydtdx)+2.*(1+A).*ux.^5.*uy.*(3.*uxdtdy+(-5).*uydtdx)+2.*(1+A).*ux.*uy.^5.*(5.*uxdtdy+(-3).*uydtdx)+8.*(1+A).*ux.*uy.*(uxdtdy+(-1).*uydtdx)+16.*(1+A).*ux.^3.*uy.^3.*(uxdtdy+(-1).*uydtdx)+(-4).*uydtdy+(-4).*A.*uydtdy+(-2).*(1+A).*ux.^6.*((-2).*uxdtdx+uydtdy)+(-26).*(1+A).*ux.^2.*uy.^2.*((-1).*uxdtdx+uydtdy)+(-2).*(1+A).*uy.^6.*((-1).*uxdtdx+2.*uydtdy)+(-2).*(1+A).*ux.^4.*((-7).*uxdtdx+4.*uydtdy)+(-2).*(1+A).*ux.^2.*((-7).*uxdtdx+5.*uydtdy)+(-2).*(1+A).*ux.^4.*uy.^2.*((-7).*uxdtdx+6.*uydtdy)+(-2).*(1+A).*ux.^2.*uy.^4.*((-6).*uxdtdx+7.*uydtdy)+(-2).*(1+A).*uy.^2.*((-5).*uxdtdx+7.*uydtdy)+(-2).*(1+A).*uy.^4.*((-4).*uxdtdx+7.*uydtdy)+(-4).*(2+A).*ux.^7.*uydt.*uydx+4.*(2+A).*ux.^3.*uy.^4.*uydt.*uydx+(-4).*ux.*(T.*uxdt+(-2).*(1+A).*uxdt.*uxdx+(3+A).*uydt.*uydx)+(-2).*ux.^5.*(T.*uxdt+(-1).*(1+A).*uxdt.*uxdx+(15+7.*A).*uydt.*uydx)+ux.^3.*((-6).*T.*uxdt+8.*(1+A).*uxdt.*uxdx+(-2).*(17+7.*A).*uydt.*uydx)+(-2).*(2+A).*ux.^6.*(1+ux.^2+uy.^2).^(1/2).*(uydt.^2+uydx.^2)+2.*(1+A).*ux.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdtdt+3.*uxdxdx+2.*uxdydy+(-5).*uydxdy)+2.*(1+A).*ux.*uy.^4.*(1+ux.^2+uy.^2).^(1/2).*(3.*uxdtdt+uxdxdx+2.*uxdydy+(-3).*uydxdy)+4.*(1+A).*ux.^3.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdtdt+uxdxdx+uxdydy+(-2).*uydxdy)+(-4).*(1+A).*ux.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdtdt+(-1).*uxdxdx+uydxdy)+(-6).*(1+A).*ux.^3.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdtdt+(-1).*uxdxdx+uydxdy)+(-2).*(1+A).*ux.^5.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdtdt+(-1).*uxdxdx+uydxdy)+4.*(2+A).*ux.^6.*uy.*(uxdt.*uydx+uydt.*(uxdx+(-1).*uydy))+(-4).*(2+A).*ux.^2.*uy.^5.*(uxdt.*uydx+uydt.*(uxdx+(-1).*uydy))+4.*(2+A).*ux.^3.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdy+uydx).*((-1).*uxdx+uydy)+4.*(2+A).*ux.^5.*uy.^2.*((-1).*uxdt.*uxdx+uxdy.*uydt+uxdt.*uydy)+(-4).*(2+A).*ux.*uy.^6.*((-1).*uxdt.*uxdx+uxdy.*uydt+uxdt.*uydy)+(-2).*ux.*uy.^2.*(5.*T.*uxdt+(-13).*uxdt.*uxdx+(-11).*A.*uxdt.*uxdx+uxdy.*uydt+(-1).*A.*uxdy.*uydt+6.*uydt.*uydx+2.*A.*uydt.*uydx+(-1).*((-1)+A).*uxdt.*uydy)+ux.^3.*uy.^2.*((-8).*T.*uxdt+8.*uxdt.*uxdx+8.*A.*uxdt.*uxdx+8.*uxdy.*uydt+8.*A.*uxdy.*uydt+(-12).*uydt.*uydx+(-4).*A.*uydt.*uydx+8.*(1+A).*uxdt.*uydy)+(-2).*ux.*uy.^4.*(3.*T.*uxdt+(-13).*uxdt.*uxdx+(-9).*A.*uxdt.*uxdx+6.*uxdy.*uydt+2.*A.*uxdy.*uydt+uydt.*uydx+A.*uydt.*uydx+2.*(3+A).*uxdt.*uydy)+(-4).*(2+A).*ux.^5.*uy.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdt.*uydt+(-1).*uxdx.*uydx+uydx.*uydy)+(-2).*(2+A).*ux.^2.*uy.^4.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdx.^2+uxdy.^2+(-1).*uydt.^2+2.*uxdy.*uydx+2.*uxdx.*uydy+(-1).*uydy.^2)+(-2).*(2+A).*ux.^4.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.^2+uxdx.^2+(-2).*uxdy.*uydx+(-1).*uydx.^2+(-2).*uxdx.*uydy+uydy.^2)+uy.^4.*(1+ux.^2+uy.^2).^(1/2).*(14.*uxdt.^2+8.*A.*uxdt.^2+uxdx.^2+A.*uxdx.^2+4.*(3+A).*uxdy.^2+(-1).*uydt.^2+(-1).*A.*uydt.^2+(-1).*uydy.^2+(-1).*A.*uydy.^2+2.*T.*((-1).*uxdx+uydy))+2.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdt.^2+2.*A.*uxdt.^2+uxdx.^2+A.*uxdx.^2+(-2).*A.*uxdy.^2+(-2).*uydt.^2+(-2).*A.*uydt.^2+2.*A.*uydx.^2+(-1).*uydy.^2+(-1).*A.*uydy.^2+2.*T.*((-1).*uxdx+uydy))+ux.^4.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.^2+A.*uxdt.^2+uxdx.^2+A.*uxdx.^2+(-14).*uydt.^2+(-8).*A.*uydt.^2+(-12).*uydx.^2+(-4).*A.*uydx.^2+(-1).*(1+A).*uydy.^2+2.*T.*((-1).*uxdx+uydy))+ux.^2.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.^2+3.*A.*uxdt.^2+6.*uxdx.^2+6.*A.*uxdx.^2+(-5).*uydt.^2+(-3).*A.*uydt.^2+(-6).*(1+A).*uydy.^2+4.*T.*((-1).*uxdx+uydy))+uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(14.*uxdt.^2+10.*A.*uxdt.^2+3.*uxdx.^2+3.*A.*uxdx.^2+(-2).*((-4)+A).*uxdy.^2+(-4).*uydt.^2+(-4).*A.*uydt.^2+2.*A.*uydx.^2+(-3).*uydy.^2+(-3).*A.*uydy.^2+6.*T.*((-1).*uxdx+uydy))+ux.^2.*(1+ux.^2+uy.^2).^(1/2).*(4.*uxdt.^2+4.*A.*uxdt.^2+3.*uxdx.^2+3.*A.*uxdx.^2+(-2).*A.*uxdy.^2+(-14).*uydt.^2+(-10).*A.*uydt.^2+(-8).*uydx.^2+2.*A.*uydx.^2+(-3).*uydy.^2+(-3).*A.*uydy.^2+6.*T.*((-1).*uxdx+uydy))+(-4).*(2+A).*ux.*uy.^5.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.*uydt+uxdy.*((-1).*uxdx+uydy))+2.*ux.^3.*uy.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.*uydt+3.*A.*uxdt.*uydt+(-2).*T.*(uxdy+(-1).*uydx)+3.*uxdx.*uydx+(-1).*A.*uxdx.*uydx+(-9).*uydx.*uydy+(-5).*A.*uydx.*uydy+(1+A).*uxdy.*(uxdx+uydy))+(-2).*ux.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.*uydt+3.*A.*uxdt.*uydt+2.*T.*(uxdy+(-1).*uydx)+uxdx.*uydx+A.*uxdx.*uydx+uydx.*uydy+A.*uydx.*uydy+(-1).*uxdy.*((9+5.*A).*uxdx+((-3)+A).*uydy))+uy.^3.*(2.*(17+7.*A).*uxdt.*uxdy+2.*uydt.*(3.*T+(-4).*(1+A).*uydy))+4.*uy.*((3+A).*uxdt.*uxdy+uydt.*(T+(-2).*(1+A).*uydy))+2.*uy.^5.*((15+7.*A).*uxdt.*uxdy+uydt.*(T+(-1).*(1+A).*uydy))+(-2).*ux.*uy.*(1+ux.^2+uy.^2).^(1/2).*(2.*T.*(uxdy+(-1).*uydx)+(-1).*uxdy.*((5+3.*A).*uxdx+(1+3.*A).*uydy)+uydx.*((1+3.*A).*uxdx+(5+3.*A).*uydy))+2.*ux.^4.*uy.*((1+A).*uxdt.*uxdy+2.*(3+A).*uxdt.*uydx+uydt.*(3.*T+2.*(3+A).*uxdx+(-1).*(13+9.*A).*uydy))+2.*ux.^2.*uy.*(2.*(3+A).*uxdt.*uxdy+(-1).*((-1)+A).*uxdt.*uydx+uydt.*(5.*T+(-1).*((-1)+A).*uxdx+(-1).*(13+11.*A).*uydy))+4.*ux.^2.*uy.^3.*((3+A).*uxdt.*uxdy+(-2).*((1+A).*uxdt.*uydx+uydt.*((-1).*T+(1+A).*(uxdx+uydy))))+(-4).*(1+A).*uy.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdxdy+uydtdt+uydydy)+(-6).*(1+A).*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdxdy+uydtdt+uydydy)+(-2).*(1+A).*uy.^5.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdxdy+uydtdt+uydydy)+(-4).*(1+A).*ux.^2.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*((-2).*uxdxdy+2.*uydtdt+uydxdx+uydydy)+(-2).*(1+A).*ux.^4.*uy.*(1+ux.^2+uy.^2).^(1/2).*((-3).*uxdxdy+3.*uydtdt+2.*uydxdx+uydydy)+(-2).*(1+A).*ux.^2.*uy.*(1+ux.^2+uy.^2).^(1/2).*((-5).*uxdxdy+5.*uydtdt+2.*uydxdx+3.*uydydy));

    PP2 = (1/4).*T.*(1+ux.^2+uy.^2).^(-3/2).*(2+ux.^2+uy.^2).^(-1).*(2.*uxdtdy+2.*A.*uxdtdy+2.*(1+A).*ux.*uxdtdx.*uy.^5+4.*(1+A).*uxdtdy.*uy.^6+(-4).*(2+A).*ux.*uxdt.*uxdy.*uy.^6+2.*(1+A).*ux.*uxdxdy.*uy.^4.*(1+ux.^2+uy.^2).^(1/2)+(-2).*(2+A).*ux.*(uxdt.^2+uxdy.^2).*uy.^5.*(1+ux.^2+uy.^2).^(1/2)+2.*(1+A).*(uxdtdt+uxdydy).*uy.^5.*(1+ux.^2+uy.^2).^(1/2)+2.*uydtdx+2.*A.*uydtdx+4.*(1+A).*ux.^6.*uydtdx+10.*(1+A).*ux.^2.*uy.^2.*(uxdtdy+uydtdx)+4.*(1+A).*uy.^2.*(2.*uxdtdy+uydtdx)+2.*(1+A).*ux.^2.*uy.^4.*(3.*uxdtdy+uydtdx)+2.*(1+A).*uy.^4.*(5.*uxdtdy+uydtdx)+4.*(1+A).*ux.^2.*(uxdtdy+2.*uydtdx)+2.*(1+A).*ux.^4.*uy.^2.*(uxdtdy+3.*uydtdx)+2.*(1+A).*ux.^4.*(uxdtdy+5.*uydtdx)+2.*(1+A).*ux.^5.*uy.*uydtdy+2.*(1+A).*ux.*uy.*(uxdtdx+uydtdy)+2.*(1+A).*ux.^3.*uy.^3.*(uxdtdx+uydtdy)+2.*(1+A).*ux.*uy.^3.*(2.*uxdtdx+uydtdy)+2.*(1+A).*ux.^3.*uy.*(uxdtdx+2.*uydtdy)+(-4).*(2+A).*ux.^6.*uy.*uydt.*uydx+(-2).*(2+A).*ux.^5.*uy.*(1+ux.^2+uy.^2).^(1/2).*(uydt.^2+uydx.^2)+2.*(1+A).*ux.^5.*(1+ux.^2+uy.^2).^(1/2).*(uydtdt+uydxdx)+2.*(1+A).*ux.*(1+ux.^2+uy.^2).^(1/2).*(uxdxdy+uydtdt+uydxdx)+2.*(1+A).*ux.^3.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(uxdxdy+uydtdt+uydxdx)+2.*(1+A).*ux.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdxdy+uydtdt+uydxdx)+2.*(1+A).*ux.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdxdy+2.*(uydtdt+uydxdx))+2.*(1+A).*ux.^4.*uy.*(1+ux.^2+uy.^2).^(1/2).*uydxdy+2.*(1+A).*uy.*(1+ux.^2+uy.^2).^(1/2).*(uxdtdt+uxdydy+uydxdy)+2.*(1+A).*ux.^2.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdtdt+uxdydy+uydxdy)+2.*(1+A).*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdtdt+2.*uxdydy+uydxdy)+2.*(1+A).*ux.^2.*uy.*(1+ux.^2+uy.^2).^(1/2).*(uxdtdt+uxdydy+2.*uydxdy)+4.*(2+A).*ux.^5.*uy.^2.*(uxdt.*uydx+uydt.*(uxdx+(-1).*uydy))+4.*(2+A).*ux.^2.*uy.^5.*((-1).*uxdt.*uxdx+uxdy.*uydt+uxdt.*uydy)+4.*(2+A).*ux.^4.*uy.^3.*((-1).*uxdt.*uxdx+uydt.*(uxdy+(-1).*uydx)+uxdt.*uydy)+ux.^4.*uy.*((-4).*(2+A).*uxdt.*uxdx+uydt.*(2.*(3+A).*uxdy+(-1).*(19+9.*A).*uydx)+2.*(3+A).*uxdt.*uydy)+uy.^5.*((-2).*T.*uxdt+(-1).*uxdt.*uxdx+A.*uxdt.*uxdx+5.*uxdy.*uydt+3.*A.*uxdy.*uydt+(5+3.*A).*uxdt.*uydy)+uy.*((-2).*T.*uxdt+(-1).*uxdt.*uxdx+A.*uxdt.*uxdx+5.*uxdy.*uydt+3.*A.*uxdy.*uydt+(-1).*uydt.*uydx+A.*uydt.*uydx+(5+3.*A).*uxdt.*uydy)+uy.^3.*((-4).*T.*uxdt+(-2).*uxdt.*uxdx+2.*A.*uxdt.*uxdx+10.*uxdy.*uydt+6.*A.*uxdy.*uydt+(-1).*uydt.*uydx+A.*uydt.*uydx+2.*(5+3.*A).*uxdt.*uydy)+ux.^2.*uy.*((-2).*T.*uxdt+(-11).*uxdt.*uxdx+(-5).*A.*uxdt.*uxdx+11.*uxdy.*uydt+5.*A.*uxdy.*uydt+(-12).*uydt.*uydx+(-4).*A.*uydt.*uydx+(11+5.*A).*uxdt.*uydy)+ux.^2.*uy.^3.*((-2).*T.*uxdt+(-19).*uxdt.*uxdx+(-9).*A.*uxdt.*uxdx+21.*uxdy.*uydt+11.*A.*uxdy.*uydt+(-9).*uydt.*uydx+(-3).*A.*uydt.*uydx+(21+11.*A).*uxdt.*uydy)+(-4).*(2+A).*ux.^4.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdt.*uydt+(-1).*uxdx.*uydx+uydx.*uydy)+ux.^4.*(1+ux.^2+uy.^2).^(1/2).*((-2).*T.*uydx+(5+3.*A).*(uxdt.*uydt+uxdx.*uydx)+((-1)+A).*uydx.*uydy)+(-2).*(2+A).*ux.^3.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.^2+uxdx.^2+uydt.^2+(-2).*uxdy.*uydx+(-2).*uxdx.*uydy+uydy.^2)+(-1).*ux.*uy.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.^2+3.*A.*uxdt.^2+2.*uxdx.^2+2.*uxdy.^2+5.*uydt.^2+3.*A.*uydt.^2+(-4).*uxdy.*uydx+2.*uydx.^2+(-4).*uxdx.*uydy+2.*uydy.^2)+(-1).*ux.^3.*uy.*(1+ux.^2+uy.^2).^(1/2).*(4.*uxdt.^2+2.*A.*uxdt.^2+2.*(2+A).*uxdx.^2+9.*uydt.^2+5.*A.*uydt.^2+(-6).*uxdy.*uydx+(-2).*A.*uxdy.*uydx+6.*uydx.^2+2.*A.*uydx.^2+(-2).*(3+A).*uxdx.*uydy+2.*uydy.^2)+(-1).*ux.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(9.*uxdt.^2+5.*A.*uxdt.^2+2.*uxdx.^2+2.*(3+A).*uxdy.^2+4.*uydt.^2+2.*A.*uydt.^2+(-2).*(3+A).*uxdy.*uydx+(-2).*(3+A).*uxdx.*uydy+2.*(2+A).*uydy.^2)+4.*(2+A).*ux.^2.*uy.^4.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.*uydt+uxdy.*((-1).*uxdx+uydy))+ux.^2.*(1+ux.^2+uy.^2).^(1/2).*(9.*uxdt.*uydt+7.*A.*uxdt.*uydt+6.*uxdx.*uydx+2.*A.*uxdx.*uydx+(-2).*T.*(uxdy+2.*uydx)+4.*A.*uydx.*uydy+(-1).*((-1)+A).*uxdy.*((-1).*uxdx+uydy))+ux.^2.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(20.*uxdt.*uydt+12.*A.*uxdt.*uydt+9.*uxdx.*uydx+3.*A.*uxdx.*uydx+(-2).*T.*(uxdy+uydx)+(-9).*uydx.*uydy+(-3).*A.*uydx.*uydy+3.*(3+A).*uxdy.*((-1).*uxdx+uydy))+(-4).*(2+A).*ux.^3.*uy.^4.*(uxdt.*uxdy+(-1).*uxdt.*uydx+uydt.*((-1).*uxdx+uydy))+(1+ux.^2+uy.^2).^(1/2).*(4.*uxdt.*uydt+4.*A.*uxdt.*uydt+uxdx.*uydx+(-1).*A.*uxdx.*uydx+(-2).*T.*(uxdy+uydx)+uydx.*uydy+3.*A.*uydx.*uydy+uxdy.*((1+3.*A).*uxdx+(-1).*((-1)+A).*uydy))+ux.^5.*((5+3.*A).*uxdt.*uydx+uydt.*((-2).*T+(5+3.*A).*uxdx+((-1)+A).*uydy))+ux.*(((-1)+A).*uxdt.*uxdy+(5+3.*A).*uxdt.*uydx+uydt.*((-2).*T+(5+3.*A).*uxdx+((-1)+A).*uydy))+uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(9.*uxdt.*uydt+7.*A.*uxdt.*uydt+uxdx.*uydx+(-1).*A.*uxdx.*uydx+(-2).*T.*(2.*uxdy+uydx)+(-1).*uydx.*uydy+A.*uydx.*uydy+2.*uxdy.*(2.*A.*uxdx+(3+A).*uydy))+uy.^4.*(1+ux.^2+uy.^2).^(1/2).*((-2).*T.*uxdy+(5+3.*A).*uxdt.*uydt+uxdy.*(((-1)+A).*uxdx+(5+3.*A).*uydy))+ux.^3.*uy.^2.*((-3).*(3+A).*uxdt.*uxdy+(21+11.*A).*uxdt.*uydx+uydt.*((-2).*T+(21+11.*A).*uxdx+(-1).*(19+9.*A).*uydy))+ux.*uy.^2.*((-4).*(3+A).*uxdt.*uxdy+(11+5.*A).*uxdt.*uydx+uydt.*((-2).*T+(-1).*(11+5.*A).*((-1).*uxdx+uydy)))+ux.^3.*(((-1)+A).*uxdt.*uxdy+2.*((5+3.*A).*uxdt.*uydx+uydt.*((-2).*T+(5+3.*A).*uxdx+((-1)+A).*uydy)))+ux.*uy.^4.*((-1).*(19+9.*A).*uxdt.*uxdy+2.*((3+A).*uxdt.*uydx+uydt.*((3+A).*uxdx+(-2).*(2+A).*uydy))));

    DP1 = (P1-P1);
    DP2 = (P2-P2);
    ep1=sqrt(sum(sum(DP1.^2)));
    ep2=sqrt(sum(sum(DP2.^2)));
%    disp(sprintf('%f %f %f %f', ep1, ep2, sqrt(sum(sum(P1.^2))), sqrt(sum(sum(P2.^2))) ));
       
    if (METRIC_ON)    
      % call the metric code
      metricz(T, ut, ux, uy, utdt, uxdt, uydt, utdx, utdy, uxdx, uxdy, uydx, uydy, ...
          utdtdt, utdtdx, utdtdy, utdxdx, utdxdy, utdydy, ...
          uxdtdt, uxdtdx, uxdtdy, uxdxdx, uxdxdy, uxdydy, ...
          uydtdt, uydtdx, uydtdy, uydxdx, uydxdy, uydydy);

%keyboard()
      % store the einstein tensor
      % res(t, :, :, :, :, :) = Gdd(:,:,:, :, :);
      
      % ************************************************
      %   check if satisfy Einstein Field Equations
      % ************************************************

      G = zeros(GRID_SIZE ,GRID_SIZE, Z_SIZE);
      % in 3+1 dim, Gdd = Rdd - (1/2R + 3) gdd == 0, so that R - (1/2R + 3) 4 == 0, so that R = -12
      % Thus G = -R - 12

      G = -gr - 12;

      % display the result
      gr = G;
      
      % integrate |G| with L2 norm;
      %  square G
      G = G.*G;

      Gsum = sqrt(sum(sum(sum(G))).*DX.*DX.*DZ)
      res(t) = Gsum;
    end

    % continue updating from EOM

    T  = T  + Dt * Tdt;
    ux = ux + Dt * uxdt;
    uy = uy + Dt * uydt;
    P1 = P1 + Dt * P1dt;
    P2 = P2 + Dt * P2dt;

    % package in res
    Nt = GRID_SIZE*GRID_SIZE;
    res(t, 0*Nt+1 : 1*Nt) = T(:);
    res(t, 1*Nt+1 : 2*Nt) = ux(:);
    res(t, 2*Nt+1 : 3*Nt) = uy(:);
    res(t, 3*Nt+1 : 4*Nt) = P1(:);
    res(t, 4*Nt+1 : 5*Nt) = P2(:);



    Tdx = dx(T);
    Tdy = dy(T);
    uxdx = dx(ux);
    uxdy = dy(ux);
    uydx = dx(uy);
    uydy = dy(uy);
    P1dx = dx(P1);
    P1dy = dy(P1);
    P2dx = dx(P2);
    P2dy = dy(P2);

% end of EOM update
    



    %if rem(t, 10) == 0
      if (SHOW_PLOT_HYDRO)
        phydro();
      end
      if (METRIC_ON && SHOW_PLOT_METRIC)
        pmetric2();
      end
      if (HORIZON_ON && SHOW_PLOT_HORIZON)
        phorizon();
      end
      if ((SHOW_PLOT_HYDRO) || (METRIC_ON && SHOW_PLOT_METRIC) || (HORIZON_ON && SHOW_PLOT_HORIZON))
        refresh();
      end
    %end
    disp(sprintf('time step: %04d, total time: %f,  step realtime: %8.4fs', t, t*Dt, toc));


  end














  function rres = fgtt(rval)
    f = (1-(T(xx,yy)./rval).^3);
    rres = -rval.^2.*f.*ut(xx,yy).*ut(xx,yy) + rval.^2.*(-1 + ut(xx,yy).*ut(xx,yy));
    delu = utdt(xx,yy)+uxdx(xx,yy)+uydy(xx,yy);
    rres = res + rval.*delu.*ut(xx,yy).*ut(xx,yy);
    rres = res -2.*rval.*ut(xx,yy).*(ut(xx,yy).*utdt(xx,yy)+ux(xx,yy).*utdx(xx,yy)+uy(xx,yy).*utdy(xx,yy));

    br = rval./T(xx,yy);
    F = -sqrt(3)./3.*atan(sqrt(3)./3.*(2.*br+1))+ 1./2.*log(1+1./br+1./(br.*br))+sqrt(3).*pi./6;
    FF = 2.*rval.^2./T(xx,yy).*F;
    rres = res + FF.*((1/2).*((-1).*utdt(xx,yy)+2.*ut(xx,yy).^2.*utdt(xx,yy)+(-1).*ut(xx,yy).^4.*utdt(xx,yy)+ut(xx,yy).*utdx(xx,yy).*ux(xx,yy)+(-1).*ut(xx,yy).^3.*utdx(xx,yy).*ux(xx,yy)+(-1).*ut(xx,yy).*ux(xx,yy).*uxdt(xx,yy)+ut(xx,yy).^3.*ux(xx,yy).*uxdt(xx,yy)+uxdx(xx,yy)+(-1).*ut(xx,yy).^2.*uxdx(xx,yy)+ux(xx,yy).^2.*uxdx(xx,yy)+ut(xx,yy).^2.*ux(xx,yy).^2.*uxdx(xx,yy)+ut(xx,yy).*utdy(xx,yy).*uy(xx,yy)+(-1).*ut(xx,yy).^3.*utdy(xx,yy).*uy(xx,yy)+ux(xx,yy).*uxdy(xx,yy).*uy(xx,yy)+ut(xx,yy).^2.*ux(xx,yy).*uxdy(xx,yy).*uy(xx,yy)+(-1).*ut(xx,yy).*uy(xx,yy).*uydt(xx,yy)+ut(xx,yy).^3.*uy(xx,yy).*uydt(xx,yy)+ux(xx,yy).*uy(xx,yy).*uydx(xx,yy)+ut(xx,yy).^2.*ux(xx,yy).*uy(xx,yy).*uydx(xx,yy)+uydy(xx,yy)+(-1).*ut(xx,yy).^2.*uydy(xx,yy)+uy(xx,yy).^2.*uydy(xx,yy)+ut(xx,yy).^2.*uy(xx,yy).^2.*uydy(xx,yy)));
    
  end

  function phorizon()
    figure(fighorizon);
    [xp,yp] = ndgrid(1:GRID_SIZE, 1:GRID_SIZE);
    size(xp)
    size(yp)
    size(horizon)
    surf(xp, yp, horizon);
    zlim( [5-4 5+4])
    caxis([5-4 5+4])
  end

  function pmetric2()
    figure(figmetric);
    plotsx = 2; 
    plotsy = 2;

    plotmag = 100;


    mdata = reshape(gr(GRID_SIZE/2, :,:), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 1);
    image(mdata*CMAP_OFFSET/10000+CMAP_OFFSET);

    subplot(plotsy, plotsx, 2);
    image(mdata*CMAP_OFFSET/100000+CMAP_OFFSET);

    subplot(plotsy, plotsx, 3);
    image(mdata*CMAP_OFFSET/1000000+CMAP_OFFSET);

    subplot(plotsy, plotsx, 4);
    image(mdata*CMAP_OFFSET/10000000+CMAP_OFFSET);

    if (VIDEO) 
      frame = getframe(gcf);
      writeVideo(vt,frame);
    end

  end
 
  function pmetric()
    figure(figmetric);
    plotsx = 3; 
    plotsy = 4;

    plotmag = 100;

    mdata = reshape(g(GRID_SIZE/2, :,:,1,1), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 1);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    title('gtt');

    mdata = reshape(g(GRID_SIZE/2, :,:,1,2), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 2);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    title('gtx');

    mdata = reshape(g(GRID_SIZE/2, :,:,1,3), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 3);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    title('gty');

    mdata = reshape(g(GRID_SIZE/2, :,:,2,2), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 4);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    title('gxx');

    mdata = reshape(g(GRID_SIZE/2, :,:,2,3), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 5);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    title('gxy');

    mdata = reshape(g(GRID_SIZE/2, :,:,3,3), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 6);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    title('gyy');

    mdata = reshape(g(GRID_SIZE/2, :,:,1,4), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 7);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    title('gtr');

    mdata = reshape(g(GRID_SIZE/2, :,:,2,4), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 8);
    image(mdata*CMAP_OFFSET+CMAP_OFFSET);
    title('gxr');

    mdata = reshape(g(GRID_SIZE/2, :,:,3,4), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 9);
    image(mdata*CMAP_OFFSET+CMAP_OFFSET);
    title('gyr');
    %if (GRID_SIZE >= 64)
    %  truesize(figmetric);
    %end
    mdata = reshape(detg(GRID_SIZE/2, :,:), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 10);
    image(mdata*CMAP_OFFSET/1000+CMAP_OFFSET);
    title('detg');

    mdata = reshape(gr(GRID_SIZE/2, :,:), GRID_SIZE, Z_SIZE);
    subplot(plotsy, plotsx, 11);
    image(mdata*CMAP_OFFSET/100000+CMAP_OFFSET);
    title('gr');

  end

  function r = phydro()
    frame_active = (mod(t-1,int16(1/Dt/SECONDS_PER_LIGHTYEAR/FPS)) == 0) || EVERY_FRAME;
    if VIDEO & ~frame_active 
      return
    end
    if VIDEO
      disp(t);
    end
    plotsx = 4;
    plotsy = 3;

    figure(fighydro);

    subplot(plotsy, plotsx, 1);
    image((T-TEMP_OFFSET)*CMAP_OFFSET +CMAP_OFFSET); 
    title('T');
    
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vt,frame);
    % end

    subplot(plotsy, plotsx, 2);
    image(ux*CMAP_OFFSET +CMAP_OFFSET);
    title('ux');
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvx,frame);
    % end

    subplot(plotsy, plotsx, 3);
    image(uy*CMAP_OFFSET +CMAP_OFFSET);
    title('uy');
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvy,frame);
    % end

    subplot(plotsy, plotsx, 4);
    image(sqrt(ux.*ux + uy.*uy)*CMAP_OFFSET +CMAP_OFFSET);
    title('sqrt(ux^2+uy^2)');
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vv,frame);
    % end

    subplot(plotsy, plotsx, 5);
    image(Tdt+CMAP_OFFSET);
    title('dtT');

    subplot(plotsy, plotsx, 6);
    image(uxdt*CMAP_OFFSET/32+CMAP_OFFSET);
    title('dtux');
    subplot(plotsy, plotsx, 7);
    image(uydt*CMAP_OFFSET/32+CMAP_OFFSET);
    title('dtuy');

    subplot(plotsy, plotsx, 9);
    image(P1*CMAP_OFFSET/8+CMAP_OFFSET);
    title('dtP1');
    subplot(plotsy, plotsx, 10);
    image(P2*CMAP_OFFSET/8+CMAP_OFFSET);
    title('dtP2');


    title(['t = ',num2str(t*Dt)]);

    if (GRID_SIZE >= 64)
      truesize(fighydro);
    end

    %if (VIDEO) & (frame_active)
    %  frame = getframe(gcf);
    %  writeVideo(vt,frame);
    %end
  end

  % [-1 0 1
  function r = dx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = - arg(GRID_SIZE,y)/2 + arg(2, y)/2;
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  -arg(GRID_SIZE-1,y)/2 + arg(1, y)/2;
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = (arg(x+1,y) - arg(x-1,y))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = - arg(x,GRID_SIZE)/2 + arg(x,2)/2;
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  -arg(x,GRID_SIZE-1)/2 + arg(x,1)/2;
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = (arg(x,y+1) - arg(x,y-1))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dxdx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = arg(GRID_SIZE,y) + arg(2, y) - 2*arg(1,y);
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  arg(GRID_SIZE-1,y) + arg(1, y) - 2*arg(GRID_SIZE,y);
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = arg(x+1,y) + arg(x-1,y) - 2*arg(x,y);
      end
    end 
    r = r / DX / DX;
  end
  function r = dydy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = arg(x,GRID_SIZE) + arg(x,2) - 2*arg(x,1);
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  arg(x,GRID_SIZE-1) + arg(x,1) - 2*arg(x,GRID_SIZE);
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = arg(x,y+1) + arg(x,y-1) - 2*arg(x,y);   
      end
    end 
    r = r / DX / DX;
  end

end


  function r = dx(arg)
    global d1;
    r = d1 * arg;
  end
  function r = dy(arg)
    global d1;
    r = arg * (d1.');
  end
  function r = dxdx(arg)
    global d2;
    r = d2 * arg;
  end
  function r = dydy(arg)
    global d2;
    r = arg * (d2.');
  end
